﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Data;
using System.Drawing;
using System.IO;
using System.Text;
using System.Windows.Forms;
using LidarReader;
using MapWindow.Map;
using MapWindow.Components;
using MapWindow.Plugins;
using MapWindow.Data;
using MapWindow.Geometries;
using Point = MapWindow.Geometries.Point;
using MapWindow.Analysis.Topology.KDTree;
using MapWindow.Drawing;
namespace HydrographicFeaturesExtract
{
    [MapWindowPlugin("HydroFeatureExtract", Author = "Ping", UniqueName = "HydroFeatureExtract", Version = "1.0")]
    public class Hydrography : Extension, IMapPlugin
    {
        #region Public
        public Raster GridFile = new Raster();
        public FeatureSet LidarFeatureSet = new FeatureSet();//correspoding shapefile
        public FeatureSet CompressedLidarPoint = new FeatureSet();//Compressed LiDAR Points
        public Raster DemRaster = new Raster();//to load the DEM file 
        public FeatureSet Outlet = new FeatureSet();//the outlet shapefile
        public int nNodes = 8000;
        public Node[] Nodes;
        public int DrawPattern = 0;
        public FeatureSet HigherNode = new FeatureSet(FeatureTypes.Point);
        public FeatureSet LineNode = new FeatureSet(FeatureTypes.Line);
        public Coordinate CurrentNode = new Coordinate();
        public KDTree _tree = new MapWindow.Analysis.Topology.KDTree.KDTree(2);
        #endregion
        private ToolStripMenuItem menuItem1;
        private ToolStripMenuItem menuItem2;
        private MapWindow.Map.IMapPluginArgs _mapWin;
        #region IMapPlugin Members
        void IMapPlugin.Initialize(IMapPluginArgs args)
        {
            _mapWin = args;
            menuItem1 = new ToolStripMenuItem("LiDAR Data Provider");
            ToolStripMenuItem menuItem2 = new ToolStripMenuItem("Load LAS Data");
            ToolStripMenuItem menuItem3 = new ToolStripMenuItem("LAS Thinning");
            ToolStripMenuItem menuItem4 = new ToolStripMenuItem("LAS2DEM");
            ToolStripMenuItem menuItem5 = new ToolStripMenuItem("LAS2SHP");
            menuItem1.DropDownItems.Add(menuItem2);
            menuItem1.DropDownItems.Add(menuItem3);
            menuItem1.DropDownItems.Add(menuItem4);
            menuItem1.DropDownItems.Add(menuItem5);
            _mapWin.MainMenu.Items.Add(menuItem1);
            //handle the menu item click event
            menuItem2.Click += new EventHandler(menuItem2_Click);
            menuItem3.Click += new EventHandler(menuItem3_Click);
            menuItem4.Click += new EventHandler(menuItem4_Click);
            menuItem5.Click += new EventHandler(menuItem5_Click);

            menuItem2 = new ToolStripMenuItem("WaterShedHydrography");
            ToolStripMenuItem menuItem6 = new ToolStripMenuItem("Load DEM");
            ToolStripMenuItem menuItem7 = new ToolStripMenuItem("Load Outlet");
            ToolStripMenuItem menuItem8 = new ToolStripMenuItem("Load LiDAR Points");
            ToolStripMenuItem menuItem9 = new ToolStripMenuItem("Get Higher Points");
            ToolStripMenuItem menuItem10 = new ToolStripMenuItem("Get Stream Flow");

            menuItem2.DropDownItems.Add(menuItem6);
            menuItem2.DropDownItems.Add(menuItem7);
            menuItem2.DropDownItems.Add(menuItem8);
            menuItem2.DropDownItems.Add(menuItem9);
            menuItem2.DropDownItems.Add(menuItem10);

            _mapWin.MainMenu.Items.Add(menuItem2);
            menuItem6.Click += new EventHandler(menuItem6_Click);
            menuItem7.Click += new EventHandler(menuItem7_Click);
            menuItem7.Click += new EventHandler(menuItem7_Click);
            menuItem8.Click += new EventHandler(menuItem8_Click);
            menuItem9.Click += new EventHandler(menuItem9_Click);
            menuItem10.Click += new EventHandler(menuItem10_Click);

            //do define a paint event which aims to draw points and lines
           // _mapWin.Map.Paint += ToDrawGraph;
        }
        #endregion
        private void menuItem2_Click(object sender, EventArgs e)//Load Lidar File(maybe from different source
        {
            ReadLasFile();
            _mapWin.Map.Layers.Add(LidarFeatureSet);
        }
        private void menuItem3_Click(object sender, EventArgs e)//to thing the file for compression
        {
            //OpenFileDialog ofd = new OpenFileDialog();
            //ofd.Filter = "Grid|*.asc";
            //if (ofd.ShowDialog() != DialogResult.OK) return;
            //GridFile.Open(ofd.FileName);
            //_mapWin.Map.Layers.Add(GridFile);
            //_mapWin.Map.Extents = GridFile.Bounds.Envelope;
            //_mapWin.Map.Layers.Add(LidarFeature);
        }
        private void menuItem4_Click(object sender, EventArgs e)//From LAS to DEM(maybe .tif)
        {
            //D8Computaion(GridFile);
            //_mapWin.Map.Extents = LidarFeature.Bounds.Envelope;
        }
        private void menuItem5_Click(object sender, EventArgs e)//From .las to .shp
        {
            OpenFileDialog ofd = new OpenFileDialog();
            ofd.Filter = "LAS|*.las";
            if (ofd.ShowDialog() != DialogResult.OK) return;
            LAS2SHP(ofd.FileName);
        }
        private void menuItem6_Click(object sender, EventArgs e)//Load DEM
        {
            DemRaster.Open("E:\\WaterShedModel\\HydroLiDAR\\Data\\iron.bgd");
            IMapRasterLayer layer = _mapWin.Map.Layers.Add(DemRaster);
            layer.Symbolizer.ShadedRelief.IsUsed = true;
            layer.Symbolizer.ShadedRelief.ElevationFactor = 1;
            layer.Symbolizer.Scheme.ApplyScheme(ColorSchemes.Summer_Mountains, DemRaster);
            layer.WriteBitmap();
            _mapWin.Map.Extents = DemRaster.Bounds.Envelope;
        }
        private void menuItem7_Click(object sender, EventArgs e)//Load Outlet
        {
            Outlet.Open("E:\\WaterShedModel\\HydroLiDAR\\Data\\ironw.shp");
            _mapWin.Map.Layers.Add(Outlet);
        }
        private void menuItem8_Click(object sender, EventArgs e)//Load LiDAR Points
        {
            //ReadXYZFile();
            //ReadLasFile();
            LidarFeatureSet.Open("E:\\WaterShedModel\\HydroLiDAR\\Data\\Lidar.shp");
            _mapWin.Map.Layers.Add(LidarFeatureSet);
            //ReadXYZFile();
            //_mapWin.Map.Layers.Add(CompressedLidarPoint);
            _mapWin.Map.Extents = LidarFeatureSet.Envelope;
        }
        private void menuItem9_Click(object sender, EventArgs e)//Get Higher Points
        {
            BuilDKDTree(LidarFeatureSet);
            ProcessNode(Nodes[1]);
            DrawPattern = 2;
            //Use a Featureset to store the polyline points of the higher points
            Point HP = new Point();
            for (int i = 0; i < Nodes.Length; i++)
            {
                if (Nodes[i].DownNode != null )
                {
                    Point[] LinePoints = new Point[2];
                    HP.X = Nodes[i].X;
                    HP.Y = Nodes[i].Y;
                    HP.Z = Nodes[i].Z;
                    LinePoints[0] = new Point(Nodes[i].X, Nodes[i].Y, Nodes[i].Z);
                    LinePoints[1] = new Point(Nodes[i].DownNode.X, Nodes[i].DownNode.Y, Nodes[i].DownNode.Z);
                    LineString ls = new LineString(LinePoints);
                    IFeature LineFeature = LineNode.AddFeature(ls);
                }
                IFeature CurrentFeature = HigherNode.AddFeature(HP);

            }
            _mapWin.Map.Layers.Add(HigherNode);
            //_mapWin.Map.Layers.Add(LineNode);
            //_mapWin.Map.Extents = LineNode.Envelope;
        }
        private void menuItem10_Click(object sender, EventArgs e)//Draw Lineation
        {
            //for (int i = 0; i < Nodes.Length; i++)
            //{
            //    if (Nodes[i].DownNode != null)
            //        DrawDelination(Nodes[i], Nodes[i].DownNode);
            //}
            _mapWin.Map.Layers.Add(LineNode);
            _mapWin.Map.Extents = HigherNode.Envelope;
        }

        private void ToDrawGraph(object sender, PaintEventArgs e)
        {
            if (Nodes == null || Nodes.Length == 0) return;
            switch (DrawPattern)
            {
                case 1://for drawing circle
                    break;
                case 2://Drawing lines
                    for (int i = 0; i < Nodes.Length; i++)
                    {
                        if (Nodes[i].DownNode != null)
                            DrawDelination(Nodes[i], Nodes[i].DownNode, e.Graphics);
                    }
                    break;
            }
        }
        private void DrawDelination(Node CurrentNode, Node DownNode, Graphics g)
        {
            System.Drawing.Point pt1 = _mapWin.Map.ProjToPixel(new Coordinate(CurrentNode.X, CurrentNode.Y));
            System.Drawing.Point pt2 = _mapWin.Map.ProjToPixel(new Coordinate(DownNode.X, DownNode.Y));
            g.FillRectangle(Brushes.Green, new Rectangle((pt1.X - 5), (pt1.Y - 5), 10, 10));
            g.DrawLine(Pens.Blue,pt1,pt2);
        }
        private void D8Computaion(Raster GridFile)
        {
            MessageBox.Show("This is the D8 Direction Algorithm");
        }
        #region IExtension Members

        void IExtension.Activate()
        {
        }

        void IExtension.Deactivate()
        {
            if (_mapWin.MainMenu.Items.Contains(menuItem1))
                _mapWin.MainMenu.Items.Remove(menuItem1);
            if (_mapWin.MainMenu.Items.Contains(menuItem2))
                _mapWin.MainMenu.Items.Remove(menuItem2);
        }

        bool IExtension.IsEnabled
        {
            get
            {
                throw new NotImplementedException();
            }
            set
            {
                throw new NotImplementedException();
            }
        }

        #endregion

        public void ReadXYZFile()//Read .XYZ file 
        {
            StreamReader sr = new StreamReader(@"E:\North_Slope_AK\Lidar_ground_AK_State_Plane_Zone4_NAD_83\5s-7e-26_grd.xyz");
            String line;
            double X, Y, Z;
            int CompressRatio = 100;
            List<Point> Listptn = new List<Point>();
            while ((line = sr.ReadLine()) != null) 
            {
                String[] strs = line.Split(' ');
                X = Convert.ToDouble(strs[0]);
                Y = Convert.ToDouble(strs[1]);
                Z = Convert.ToDouble(strs[2]);
                Point pntLidar = new Point();
                pntLidar.X = X;
                pntLidar.Y = Y;
                pntLidar.Z = Z;
                Listptn.Add(pntLidar);
                MapWindow.Data.IFeature f = LidarFeatureSet.AddFeature(pntLidar);//Add the feature to the FeatureSet

                //The following is do some compressing
                if (Listptn.Count % CompressRatio == 0 && Listptn.Count != 0)
                {
                    Point CompressPnt = new Point();
                    double XCom, YCom, ZCom;
                    XCom = YCom = ZCom = 0;
                    for (int i = 0; i < Listptn.Count ; i++)
                    {
                        XCom += Listptn[i].X;
                        YCom += Listptn[i].Y;
                        ZCom += Listptn[i].Z;
                    }
                    CompressPnt.X = XCom / (double)Listptn.Count;
                    CompressPnt.Y = YCom / (double) Listptn.Count;
                    CompressPnt.Z = ZCom / (double)Listptn.Count;
                    MapWindow.Data.IFeature Pf = CompressedLidarPoint.AddFeature(CompressPnt);
                    Listptn.Clear();
                    //Listptn.Count = 0;
                }
            }
            MessageBox.Show("ok ");
            //Extract each 30 points
        }

        public void LAS2SHP(string FileName)
        {
            Stream stream = new FileStream(FileName, FileMode.Open, FileAccess.Read, FileShare.Read);
            BinaryReader br = new BinaryReader(stream);
            LasHeader hdr = new LasHeader();
            hdr.Read(br);
            br.BaseStream.Seek(hdr.PointDataOffset, SeekOrigin.Begin);
            //Loop through and read each of the points in  the LAS file
            #region //using the table to store the information read from the LiDAR Point
            LidarFeatureSet.DataTable.Columns.Add(new MapWindow.Data.Field("Id", MapWindow.Data.FieldDataTypes.String));
            LidarFeatureSet.DataTable.Columns.Add(new MapWindow.Data.Field("Blue", MapWindow.Data.FieldDataTypes.String));
            LidarFeatureSet.DataTable.Columns.Add(new MapWindow.Data.Field("EdgeFlight", MapWindow.Data.FieldDataTypes.String));
            LidarFeatureSet.DataTable.Columns.Add(new MapWindow.Data.Field("GPSTime", MapWindow.Data.FieldDataTypes.Double));
            LidarFeatureSet.DataTable.Columns.Add(new MapWindow.Data.Field("Green", MapWindow.Data.FieldDataTypes.String));
            LidarFeatureSet.DataTable.Columns.Add(new MapWindow.Data.Field("Intensity", MapWindow.Data.FieldDataTypes.String));
            LidarFeatureSet.DataTable.Columns.Add(new MapWindow.Data.Field("PtClass", MapWindow.Data.FieldDataTypes.String));
            LidarFeatureSet.DataTable.Columns.Add(new MapWindow.Data.Field("Red", MapWindow.Data.FieldDataTypes.String));
            LidarFeatureSet.DataTable.Columns.Add(new MapWindow.Data.Field("ReturnNum", MapWindow.Data.FieldDataTypes.String));
            LidarFeatureSet.DataTable.Columns.Add(new MapWindow.Data.Field("ScanAngle", MapWindow.Data.FieldDataTypes.String));
            LidarFeatureSet.DataTable.Columns.Add(new MapWindow.Data.Field("ScanDir", MapWindow.Data.FieldDataTypes.String));
            LidarFeatureSet.DataTable.Columns.Add(new MapWindow.Data.Field("SourceID", MapWindow.Data.FieldDataTypes.String));
            LidarFeatureSet.DataTable.Columns.Add(new MapWindow.Data.Field("TotalReturns", MapWindow.Data.FieldDataTypes.String));
            LidarFeatureSet.DataTable.Columns.Add(new MapWindow.Data.Field("UserData", MapWindow.Data.FieldDataTypes.String));
            #endregion
            LidarFeatureSet.FillAttributes();
            for (int i = 0; i < hdr.Count; i++)
            {
                LidarPoint lidpnt = new LidarPoint();
                lidpnt.ReadPoint(br, hdr);
                //if (lidpnt.PtClass == 2)//for the ground return
                Point pntLidar = new Point();

                pntLidar.X = lidpnt.X;
                pntLidar.Y = lidpnt.Y;
                pntLidar.Z = lidpnt.Z;
                MapWindow.Data.IFeature f = LidarFeatureSet.AddFeature(pntLidar);//Add the feature to the FeatureSet
                //****************************************************//
                #region
                //input the attribute to the table
                f.DataRow["Id"] = i.ToString();
                f.DataRow["Blue"] = lidpnt.Blue.ToString();
                f.DataRow["EdgeFlight"] = lidpnt.EdgeFlight.ToString();
                f.DataRow["GPSTime"] = lidpnt.GPSTime;
                f.DataRow["Green"] = lidpnt.Green.ToString();
                f.DataRow["Intensity"] = lidpnt.Intensity.ToString();
                f.DataRow["PtClass"] = lidpnt.PtClass.ToString();
                f.DataRow["Red"] = lidpnt.Red.ToString();
                f.DataRow["ReturnNum"] = lidpnt.ReturnNum.ToString();
                f.DataRow["ScanAngle"] = lidpnt.ScanAngle.ToString();
                f.DataRow["ScanDir"] = lidpnt.ScanDir.ToString();
                f.DataRow["SourceID"] = lidpnt.SourceID.ToString();
                f.DataRow["TotalReturns"] = lidpnt.TotalReturns.ToString();
                f.DataRow["UserData"] = lidpnt.UserData.ToString();
                #endregion
                //LidarFeatureSet.DataTable
            }
            SaveFileDialog sfd = new SaveFileDialog();
            sfd.Filter = "SHP|*.shp";
            sfd.Title = " Save to a shapefile";
            sfd.ShowDialog();
            string ShpFileName = sfd.FileName;
            LidarFeatureSet.SaveAs(ShpFileName, false);
        }
        public void ReadLasFile()
        {
            OpenFileDialog ofd = new OpenFileDialog();
            if (ofd.ShowDialog() != DialogResult.OK) return;
            Stream stream = new FileStream(ofd.FileName, FileMode.Open, FileAccess.Read, FileShare.Read);
            BinaryReader br = new BinaryReader(stream);
            LasHeader hdr = new LasHeader();
            hdr.Read(br);
            //Jump to the start of the point data
            br.BaseStream.Seek(hdr.PointDataOffset, SeekOrigin.Begin);
            int nNodes = (int)hdr.Count;
            Nodes = new Node[nNodes];
            //List<Node> Nodes = new List<Node>();     
            for (int i = 0; i < hdr.Count; i++)
            {
                LidarPoint lidpnt = new LidarPoint();
                lidpnt.ReadPoint(br, hdr);
                //if (lidpnt.PtClass == 2)//for the ground return
                Point pntLidar = new Point();

                pntLidar.X = lidpnt.X;
                pntLidar.Y = lidpnt.Y;
                pntLidar.Z = lidpnt.Z;
                MapWindow.Data.IFeature f = LidarFeatureSet.AddFeature(pntLidar);//Add the feature to the FeatureSet

                //*********To build up the tree***********************//
                Nodes[i] = new Node();
                Nodes[i].X = pntLidar.X;
                Nodes[i].Y = pntLidar.Y;
                Nodes[i].Z = pntLidar.Z;
                double[] dbl = new double[2];
                dbl[0] = Nodes[i].X;
                dbl[1] = Nodes[i].Y;
                _tree.Insert(dbl, Nodes[i]);
            }
            //ProcessNode(Nodes[100]);
          }

        private void BuilDKDTree(FeatureSet fs)
        {
            int nNodes = fs.Features.Count;
            Nodes = new Node[nNodes];
            for(int i=0;i< nNodes ;i++)
            {
                Nodes[i] = new Node();
                Nodes[i].X = fs.Features.ElementAt(i).Coordinates[0].X ;
                Nodes[i].Y = fs.Features.ElementAt(i).Coordinates[0].Y ;
                Nodes[i].Z = fs.Features.ElementAt(i).Coordinates[0].Z ;
                double[] dbl = new double[2];
                dbl[0] = Nodes[i].X;
                dbl[1] = Nodes[i].Y;
                _tree.Insert(dbl, Nodes[i]);           
            }
        }
        private void ProcessNode(Node CurrentNode)
        {
            //compute distance to every node
            //Dim d(nNodes - 1) As Double, keys(nNodes - 1) As Integer
            //double d(nNode-1);
            double [] d = new double[nNodes -1];
            int [] keys = new int[nNodes -1];
            double slope;
            double [] Location = new double[2];
            Location[0] = CurrentNode.X;
            Location[1] = CurrentNode.Y;
            //Dim Neighbors() As Object
            int K = 20;
            object[] Neighbors = new object[K];
            Neighbors = _tree.Nearest(Location, K);
            //'check nearest K nodes to see if any are higher in elevation than the current node.
            for (int i = 0; i < K; i++)
            {
                Node Neighbor = new Node();
                Neighbor = (Node)(Neighbors[i]);
                if (Neighbor.Z > CurrentNode.Z)//compute slope to this uphill node and see if it is steeper than any existing slope stored at that node
                {
                    double oldSlope = Neighbor.Slope();
                    slope = Neighbor.GetSlope(CurrentNode);           
                    if (slope > oldSlope)
                    {
                        Neighbor.DownNode = CurrentNode;
                    }
                    //now recursively do this uphill node... 
                    if (Neighbor.Visited == false) ProcessNode(Neighbor);
                }
            }
            CurrentNode.Visited = true;
            Point HigherPoint = new Point();
            HigherPoint.X = CurrentNode.X;
            HigherPoint.Y = CurrentNode.Y;
            HigherPoint.Z = CurrentNode.Z;
            MapWindow.Data.IFeature f = HigherNode.AddFeature(HigherPoint);//Add the feature to the FeatureSet
            //g.DrawEllipse(Pens.Blue, (int)(CurrentNode.X - 1), (int)(CurrentNode.Y - 1), (int)(CurrentNode.X + 1), (int)(CurrentNode.Y + 1));
            //g.DrawEllipse(Pens.Blue, RECT);
            //draw this node on the map
            //AxMap1.DrawCircle(CurrentNode.X, CurrentNode.Y, 4, cBlue, False)
            //MapWindow.PluginInterfaces.IDraw DrawCircle = new MapWindow.PluginInterfaces.IDraw();
            //DrawCircle.DrawCircle(CurrentNode.X, CurrentNode.Y, 2, Pens.Blue, false);
            //System.Drawing.Graphics g = new System.Drawing.Graphics();
            //Graphics g = map1.CreateGraphics();
            //_mapWin.Map.AddLayer();
        }
    }
    public class Node
    {
        public double X,Y,Z;
        public Node DownNode;
        public bool Visited;
        public double Slope()
        {
            if (DownNode==null)
                return 0;
            else return GetSlope(DownNode);
        }
        public double GetSlope(Node DownNode)
        {
            double dx, dy, dist;;
            dx= X - DownNode.X;
            dy= Y - DownNode.Y;
            dist =Math.Sqrt(dx * dx + dy * dy);
            return (Z - DownNode.Z) / dist;
        }
    }
    public class Zone
    {
        public int ID,Column,Row;
        public int ContaintedNodes() { return 0; }
    }
}
